;-------------------------------------------------------------LICENSE--------------------------------------------------------------;
;                                                                                                                                  ;
;The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    ;
;and Message Passing Interface (MPI) parallelization.                                                                              ;
;                                                                                                                                  ;
;Copyright (C) 2012                                                                                                                ;
;Ronglin Jiang                                                                                                                     ;
;rljiang@ssc.net.cn                                                                                                                ;
;585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             ;
;                                                                                                                                  ;
;This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         ;
;as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             ;
;                                                                                                                                  ;
;This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    ;
;MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        ;
;                                                                                                                                  ;
;You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     ;
;Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   ;
;                                                                                                                                  ;
;-------------------------------------------------------------LICENSE--------------------------------------------------------------;

pro vector_potential_2d, bx, by, x, y, a, DIR = dir, POS_START = pos_start

   if n_elements(bx) ne n_elements(by) and n_elements(bx) ne n_elements(x) * n_elements(y)   $
      then message, ' Error: bx and by must have the same size.',   $ 
                    ' x and y must have the same one dimensional size of bx and by in x and y direction.'

   nx = n_elements(x)
   ny = n_elements(y)

   if n_elements(dir) eq 0 then dir = 1
   if dir eq 1 then begin
      if not keyword_set(pos_start) then pos_start = nx / 2
   endif
   if dir eq 2 then begin
      if not keyword_set(pos_start) then pos_start = ny / 2
   endif

   bx_h = bx[*, 0 : ny - 2]
   by_h = by[0 : nx - 2, *]
   bx_h = (bx[*, 0 : ny - 2] + bx[*, 1 : ny - 1]) / 2.0
   by_h = (by[0 : nx - 2, *] + by[1 : nx - 1, *]) / 2.0

   a = bx
   a[*] = 0.0

   if dir eq 1 then begin
      for j = 0, ny - 2 do begin
         jp1 = j + 1
         a[pos_start, jp1] = a[pos_start, j] + (y[jp1] - y[j]) * bx_h[pos_start, j]
      endfor
      for j = 0, ny - 1 do begin
         for i = pos_start, nx - 2 do begin
            ip1 = i + 1
            a[ip1, j] = a[i, j] - (x[ip1] - x[i]) * by_h[i, j]
         endfor
         for i = pos_start, 1, -1 do begin
            im1 = i - 1
            a[im1, j] = a[i, j] + (x[i] - x[im1]) * by_h[i, j]
         endfor
      endfor
   endif

   if dir eq 2 then begin
      for i = 0, nx - 2 do begin
         ip1 = i + 1
         a[ip1, pos_start] = a[i, pos_start] + (x[ip1] - x[i]) * by_h[i, pos_start]
      endfor
      for i = 0, nx - 1 do begin
         for j = pos_start, ny - 2 do begin
            jp1 = j + 1
            a[i, jp1] = a[i, j] - (y[jp1] - y[j]) * bx_h[i, j]
         endfor
         for j = pos_start, 1, -1 do begin
            jm1 = j - 1
            a[i, jm1] = a[i, j] + (y[j] - y[jm1]) * bx_h[i, j]
         endfor
      endfor
   endif

   if dir eq 3 then begin

      a1 = a
      if not keyword_set(pos_start) then pos_start = nx / 2

      for j = 0, ny - 2 do begin
         jp1 = j + 1
         a1[pos_start, jp1] = a1[pos_start, j] + (y[jp1] - y[j]) * bx_h[pos_start, j]
      endfor
      for j = 0, ny - 1 do begin
         for i = pos_start, nx - 2 do begin
            ip1 = i + 1
            a1[ip1, j] = a1[i, j] - (x[ip1] - x[i]) * by_h[i, j]
         endfor
         for i = pos_start, 1, -1 do begin
            im1 = i - 1
            a1[im1, j] = a1[i, j] + (x[i] - x[im1]) * by_h[i, j]
         endfor
      endfor

      a2 = a1
      if not keyword_set(pos_start) then pos_start = ny / 2

      for i = 0, nx - 2 do begin
         ip1 = i + 1
         a2[ip1, pos_start] = a2[i, pos_start] + (x[ip1] - x[i]) * by_h[i, pos_start]
      endfor
      for i = 0, nx - 1 do begin
         for j = pos_start, ny - 2 do begin
            jp1 = j + 1
            a2[i, jp1] = a2[i, j] - (y[jp1] - y[j]) * bx_h[i, j]
         endfor
         for j = pos_start, 1, -1 do begin
            jm1 = j - 1
            a2[i, jm1] = a2[i, j] + (y[j] - y[jm1]) * bx_h[i, j]
         endfor
      endfor

      a = (a1 + a2) / 2.0
;a = a2
;a = a1

   endif

end
